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The initial trial wave function used in a simple ground-state projection 
method, the power method, is systematically improved by using Lanczos algo- 
rithm. Much faster convergence to the ground state achieved by using these 
wave functions significantly reduces the effect of the fermion sign problem. The 
energy, spin and charge correlation functions are calculated for the ground states 
of the two-dimensional t — J model. Results for an 8 x 8 cluster with 42 and 26 
electrons are presented. The density correlation function for the t — J model 
at small J shows a surprisingly good agreement with that of a system of non- 
interacting hard-core bosons. 
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Recently we have studied ground-state properties of the t — J model in one 
and two dimensions 2 by using a simplified Green function Monte Carlo (GFMC) 
method 3 - the power method. In this method the ground-state wave function of 
a Hamiltonian H is obtained by applying large powers of the operator W — H 
to a trial wave function, where W is a constant. In fermionic systems when the 
power becomes large same configurations with opposite signs will be generated 
if a Monte Carlo (MC) algorithm is used. It causes very large error bars in 
numerical values. This is the famous sign problem 3 ' 4 occurred in MC simula- 
tions of fermionic systems. In one dimension the phase of the wave function 
can be fixed to rid of the sign problem, power method is very successful for all 
possible electronic densities 1 . In two dimensions only at low electronic density 
the sign problem is not severe and the converged ground state is obtained 2 . At 
high density the sign problem makes the power method ineffective to study this 
interesting region for high-temperature superconductors. 

The freedom to choose the trial wave function is one of the special prop- 
erties of the ground-state projection method. A trial function chosen inap- 
propriately would require a lot of computer time to converge to ground state. 
Sometimes the sign problem makes the convergence impossible. It is impera- 
tive to have a good trial function to reduce the number of negative terms which 
increases with the power. 

In the last several years variational MC method has been widely used 
to study the t — J model 5 ' 6 7 . Several innovative wave functions have been 
proposed for the ground state. Some of them tested by the power method are 
not as close to the ground state as one would have anticipated. There were 
few methods that we can use to systematically improve the trial wave function. 
Recently Heeb and Rice 8 proposed to use Lanczos 9 iteration to obtain better 
wave functions. The effectiveness of the method is demonstrated by studying 
the two-dimensional antiferromagnetic Heisenberg model. A few years earlier, 
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Caffarel et al have also used a variation of the Lanczos algorithm to study 
LiH molecule. 

Although Lanczos method 9 is best known in searching for wave functions of 
small clusters, the method itself is quite general. Starting with a wave function 
|0o)> we can generate a tri-diagonal matrix by using the recurrence relation 9 

H\(f)o) = a |0 O > + &i |0i> 

H \4> n ) = a n \4> n ) + b n \4>n-l) + K+l \4>n+l) (i) 

where n =1,2, ..etc. The matrix elements, a n and b n , are related to the moments 
of the Hamiltonian. For example, ao = (<j>o \H\ 0o) an d bi = a/ (4>o \ (H — ao) 2 | (fio) 
When n increases, the lowest eigenvalue of the tri-diagonal matrix approaches 
the ground-state energy. And the eigenfunction of this lowest eigenvalue gets 
closer to the ground-state wave function. It is straight forward to show that in 
succesive iteration the eigenstates have the form 

|*i) = |0o) + Ci^|0o>, (2) 

and 

|* 2 ) = |0o) + C[±H |^o) + C^H 2 |0 O > , (3) 

... etc. These functions form the basis in Krylov subspace 11 . The C's are 
calculated from the the matrix elements, a n and 6 n , by diagonalizing the matrix. 

Heeb and Rice 8 propose to calculate the matrix elements, a n and 6 n , by 
using the Monte Carlo technique. The C's are then determined. However, in 
this method the values of the matrix elements must be calculated very accu- 
rately. A small error will produce large uncertainty in the eigenvalues and in 
C's. Here we choose an alternative. We treat C's as the variational parame- 
ters. The wave function with the optimal energy is the eigenfunction with the 
lowest eigenvalue. This is more efficient and sometimes more accurate than 
diagonalizing the matrix 8 ' 10 . 
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The result of this variational Lanczos algorithm is that we have a sequence 
of wave functions, \4>o) , |\l>i) , |*2)r- etc., with lower and lower energy. Besides 
the statistical fluctuation associated with the MC technique, the same result as 
the Lanczos method will be obtained. The fact that this method does not need 
very large memory space to store all the configurations as in the usual Lanczos 
method is one of its biggest advantages. But there is a practical difficulty 
with this approach of getting the ground state. Each time the Hamiltonian 
H is applied to a particular configuration the number of new configurations 
generated is of order of, N, the size of the cluster. It is impractical to do any 
calculation with \^ n ) for n > 3 for a cluster of 64 sites or greater. A more 
efficient way to obtain the ground state is to use or |^ 2 ) as the trial wave 
functions in the power method. We shall refer to this as the power-Lanczos 
(PL) method. If the starting trial function before the power method is applied 
is |^ n ) we shall call it PLn. PLO is the same as the usual power method. For 
the reason discussed above we shall only consider PL1 and PL2 in this paper. 

Once the optimal wave functions and j^) are determined, we can 
proceed to calculate quantities such as ^^^fj^y^ , where p is the power. 
It is sufficient to choose the constant W to be zero in the t — J model. The 
procedure to carry out this part is the same as the power method 1 . 

We use several different forms of |(/>o) to study the t — J Hamiltonian. The 
familiar Gutzwiller wave function (GWF) 12 is just the wave function for an 
ideal Fermi gas excluding configurations with doubly occupied sites. Another 
function proposed by Hellberg and Mele 6 and used by Valenti and Gros 7 in 
2D was shown to be close to the ground state at low density 2 . This function, 
which we shall call HMVG, is basically of the same form as GWF, i.e. a Slater 
determinant for up-spin electrons and one for down-spin electrons. In addition 
to these two determinants, it has a long range correlation part between all 
the particles, — Yj\ u ( while for nearest-neighbor particles we choose 
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v = 0). Besides these two functions we also use the projected BCS state or 
the resonating-valence-bond state 5 ' 13 with either s-wave or d-wave symmetry 
for the gap order parameter. 

The energy, E = (i7 2p+1 ) / (-£/" 2p ), as a function of power p is plotted in 
Fig.l for 10 particles in a 4x4 lattice. Here we consider J = 2t and GWF is 
chosen to be the initial trial function \<po). The open triangles are the result of 
Lanczos algorithm for different orders of iteration. These results are obtained 
exactly using the usual Lanczos method described briefly in equation 1. The 
variational energy of GWF is about 5% above the ground state energy. This 
difference is reduced to about 0.3% by using the second order wave function. 
The solid circles, squares and triangles are the results of PLO, PL1 and PL2 by 
using |0o) 7 1^1)7 an d ^2) respectively. For of equation (2), we choose C\ 
to be 0.8. We have C[ = 1.72 and C' 2 = 0.72 in |* 2 ) of equation (3). Clearly, 
when the power becomes large enough, all these three algorithms would produce 
ground-state energy. For comparison, we also calculated the energy exactly 
without using the MC technique in PLO, PL1 and PL2. They are shown by 
the dashed lines. The excellent agreement between exact and MC calculations 
reaffirms the stability of MC technique. 

In Fig. 1 we note the relatively large error bars at powers greater than 
4. This is mainly due to the fermion sign problem. The effect of this sign is 
studied by calculating the ratio of contributions from the negative terms and 
contributions from the positive terms in the quantities ((—H) m ). In the inset 
of figure 1 this ratio is plotted as a function of m. This ratio is about the same 
for different PL's. At power equal to 6 or m = 13, the negative terms are as 
large as seventy percent of the positive terms. It is very time consuming to get 
good statistics. 

The data for PLO and PL1 are obtained by averaging ten to twenty in- 
dependent groups. Each group usually consists of one to two thousand start- 
ing configurations. Each starting configuration would produce several hundred 
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terms in the evaluation of powers of H. At large power we usually need to have 
more runs to reduce the fluctuations caused by the sign problem. The calcula- 
tions are all done in workstations like HP720. For PLO, the longest calculation 
is about six hours for this 4x4 cluster. It takes about four to six times longer 
for PL1. If we had used the same number of runs for PL2 we would probably 
increase the CPU time by another factor of four to six. Instead, we reduced 
the number of runs in PL2, therefore we have a somewhat larger error at power 
equal to four. The amount of computer time quoted above is only for a small 
cluster of size 4x4. We estimate the time needed for an 8x8 cluster with same 
electronic density is about ten times longer. Clearly, it is quite impractical to 
carry out a calculation for large powers in PL2. Hence, below we shall present 
mostly results of PLO and PL1 for 8x8 clusters and may be a few results of 
PL2 without powers. 

In figure 2(a) and 3(a) energy as a function of power is plotted for 26 and 
42 electrons repectively, in an 8 x 8 lattice for J = O.lt. VGHM function with 
v = 0.04 is used as |0 O ) and its variational energy is about 3% above the ground 
state energy for (n) = ||, but more than 5% for (n) = ||. The situation is 
improved substantially in PL1 when of equation (2) is used. C\ = 1.33 in 
2(a) and 1.49 in 3(a). In figure 2(a) the large error bars at power equal to 6 
and 8 makes it difficult to determine the exact ground state energy. The effect 
of the negative sign is larger when the density is increased from (n) = || to ||. 
Neverthless, it seems quite reasonable to conclude that the ground-state energy 
should be within half a percent of the variational energy of j^) which are 
represented by the solid triangles. For PL2 in Fig. 2, C[ = 2.66 and C 2 = 1-77, 
and C[ = 2.95 and C 2 = 2.18 in Fig.3. 

Besides the energy we also calculate the equal time correlation functions, in 
particular, the spin and density structure factors, S(k) and N(k) respectively. 
These structure factors are plotted along r-X-M-T direction in the Brillouin 
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zone in figures 2(b) and 2(c) for (n) = g| and in figures 3(b) and 3(c) for 
(n) = || . Open circles represent the variational result of VGHM jy=0 .04 5 and 
open squares are for PL1 without power. Open triangles are results of PL1 
with power equal to 6 for (n) = || and power equal to 5 for (n) = ||. The 
solid lines connecting triangles are guides for the eyes. We note that the results 
changed markedly between the initial variational wave function and the first 
order Lanczos wave function. The situation seems to get worse when the density 
(n) increases, even though VGHM wave function still has the best variational 
energy at J = O.lt. This points out a possible deficiency in using the trial 
wave function VGHM to understand ground states of the t — J model at high 
electronic density. For comparison, we also show the results of GWF as dotted 
lines. GWF clearly does not reflect the correlation of the ground state. So far 
we have not yet found a wave function that would have energy within 5% of 
the ground state. 

Recently we have shown that many results of the t — J model at low elec- 
tronic density are qualitatively consistent with the prediction of the Tomonaga- 
Luttinger liquid 14 in one dimension. The cusps or peaks at k = 2kp in S(k) 
are enhanced over the variational results of VGHM. N(k) has a maximum at 
k = (it, it). But we cannot 15 identify in iV(k) the characteristic wave vector 
2kf, F associated with spinless fermions (SF) as claimed by Putikka et al. 16 us- 
ing the high temperature expansion. Here we try to understand iV(k) from a 
different point of view. 

One way to treat the constraint of no double occupancy in the t — J model 
is to write the fermion operator as a product of a hard-core boson and a fermion 
operator. While the fermion operator represents spin degree of freedom, the bo- 
son is for charge degree of freedom. This is the so called slave-boson approach 17 . 
If the separation of charge and spin indeed occurs it will be most apparent in 
the limit of vanishing J where the dynamics is controlled by the charge hopping 
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or the motion of the hard-core bosons. We may expect the charge correlation 
to be simialr to that of a system of non-interacting hard-core bosons. In one 
dimension hard-core bosons and spinless fermions are equivalent, but they are 
not in two dimensions. A careful examination of the correlation function of the 
hard-core bosons in a 2D lattice is needed. 

We have calculated the ground-state correlation function of a system of 
non-interacting hard-core bosons by using the power method. The trial wave 
function is of the form of Jastrow type described in Ref. 18. Details of this 
calculation will be presented elsewhere. Results of density correlation are rep- 
resented by the solid circles in figure 2(c) and 3(c). They almost lie exactly on 
top of the triangles representing the result of PL1, except 19 at very small k. 
A similar result 20 has been found for the infinite-U Hubbard model for small 
clusters. On the other hand, N(k) of SF as shown by the dashed lines in figure 
2(c) and 3(c) is not as close to the result oft — J model. The fact that hard-core 
bosons have almost the same density-density correlation as the charges in the 
t — J model does not by itself prove the separation of spin and charge. But this 
and other evidences 2 ' 15 ' 16 make the idea 21 of separation of charge and spin in 
the t — J model much more plausible. 

In summary, we have presented a new modified power method using a 
systematically improved trial wave function obtained by the Lanczos method. 
Even the wave function obtained by a first order iteration greatly improves the 
rate of convergence to the ground state. This faster convergence significantly 
reduces the effect of the fermion sign problem that has plagued fermion MC 
calculations so far. For the first time very accurate ground-state results are 
obtained for electronic density as high as 65% in an 8 x 8 cluster. A surprising 
result has been found. The density-density correlation obtained at small J is 
very close to that of a system of non-interacting hard-core bosons. 

The power-Lanczos method presented above is in principle an exact ap- 
proach to obtain the ground state. Unlike the fixed-node method 22 , the results 

8 



are not overwhelmingly influenced by the initial choice of trial function. The 
numerical algorithm we have used is the simplest among many complicated 
GFMC methods. It is already sufficient to get very accurate results. More 
sophisticated approaches, such as using guiding function 3 ' 22 , would be explored 
in the future. Calculations for much larger clusters are feasible now. 
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the National Science Council of Republic of China, Grant Nos. NSC83-0511- 
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Figure Captions: 

Fig. 1 Energy as a function of power for 10 electrons in a 4 x 4 cluster. GWF is 
the trial function used for J = 2t. The solid circles, squares and triangles 
represent results for PLO, PL1 and PL2, respectively. The dashed lines 
represent exact results without using Monte Carlo technique. Open trian- 
gles are the exact results obtained from each order of Lanczos iteration. 
In the inset ratio of contributions from negative terms and contributions 
from positive terms as a function of the power m in ((—H) m ). 

Fig. 2 (a)Energy as a function of power for (n) = || and J = O.lt calculated 
using PLO, PL1 and PL2 algorithms. Ci=1.33 in C[ = 2.66 and 

C' 2 = 1-77 in 1^2)- (b) spin structure factor S(k) and (c) density struc- 
ture factor N(k) in the k space along T-X-M-r directions. Empty circles 
represent variational results using VGHM function with v = 0.04. Open 
squares are results of PL1 without power. Open triangles are PL1 results 
at power equal to 6. Solid circles represent results of non-interacting hard- 
core bosons. Dotted line is the variational results of GWF. Results of SF 
are represented by the dashed line. 

Fig. 3 (a)Energy as a function of power for (n) = || and J = O.lt calculated 
using PLO, PL1 and PL2 algorithms. Ci=1.49 in C[ = 2.95 and 

C 2 = 2.18 in 1*2)- (b) spin structure factor S(k) and (c) density struc- 
ture factor N(k) in the k space along T-X-M-r directions. Empty circles 
represent variational results using VGHM function with v = 0.04. Open 
squares are results of PL1 without power. Open triangles are PL1 results 
at power equal to 5. Solid circles represent results of non-interacting hard- 
core bosons. Dotted line is the variational results of GWF. Results of SF 
are represented by the dashed line. 
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